function Features = produceDeepFeatures(NetFile, WeightsFile, layer, MinBatchSize, ImgIndex, hImgRead, hImgProc)

%%
mycaffe.reset_all();

if ~mycaffe.useGPU
    warning('Using CPU mode!');
    caffe.set_mode_cpu();
else
    caffe.set_mode_gpu();
    gpu_id = 0; % use the second gpu
    caffe.set_device(gpu_id);
end

net = caffe.Net(NetFile, WeightsFile, 'test');
% ImgMean = caffe.read_mean(MeanFile);
InputDim = net.blobs(layer.data).shape;

%%
if isvector(ImgIndex) && ~isstruct(ImgIndex)
    bGivenData = false;
    SampleNum = length(ImgIndex);
else
    bGivenData = true;
    if isstruct(ImgIndex)
        ImgIndex = ImgIndex.X;
    end
    X = hImgProc(ImgIndex, InputDim);
    [~,~,~,SampleNum] = size(X);
end

%%
MinBatchSize = min([MinBatchSize SampleNum]);
BatchSize = InputDim(end);
if BatchSize < MinBatchSize
    BatchSize = MinBatchSize;
    InputDim(end) = MinBatchSize;
    net.blobs(layer.data).reshape(InputDim);
end

%%
FeatureShape = net.blobs(layer.feature).shape;
% Features = zeros(LayerShape(1), SampleNum);
FeatureDim = prod(FeatureShape(1:end-1));
Features = zeros(FeatureDim, SampleNum);

%%
BatchNum = floor(SampleNum/BatchSize);
remainder = SampleNum - BatchNum * BatchSize;
if remainder > 0
    BatchNum = BatchNum + 1;
end

%%
pp = 0;
for n = 1 : BatchNum
    fprintf('%d/%d\n', n, BatchNum);
    if n == BatchNum && remainder > 0
        BatchSize = remainder;
        InputDim(end) = BatchSize;
        net.blobs(layer.data).reshape(InputDim);
    end
    
    pp = pp(end)+1 : pp(end)+BatchSize;
    
    if bGivenData
        J = X(:,:,:,pp);
    else
        X = hImgRead(ImgIndex(pp));
        J = hImgProc(X, InputDim);
    end
    
    net.forward({J});
    f = net.blobs(layer.feature).get_data();
    Features(:, pp) = reshape(f, [FeatureDim BatchSize]);
    
    % layer_conv52 = net.blob_vec(net.name2blob_index('pool5'));
    % conv52 = layer_conv52.get_data();
    % sum(conv52(:)>0) /320/100
end

